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Abstract 

We present Bill2d, a modern and efficient C++ package for classical simulations of two-dimensional Hamil¬ 
tonian systems. Bill2d can be used for various billiard and diffusion problems with one or more charged 
particles with interactions, different external potentials, an external magnetic field, periodic and open bound¬ 
aries, etc. The software package can also calculate many key quantities in complex systems such as Poincare 
sections, survival probabilities, and diffusion coefficients. While aiming at a large class of applicable systems, 
the code also strives for ease-of-use, efficiency, and modularity for the implementation of additional features. 
The package comes along with a user guide, a developer’s manual, and a documentation of the application 
program interface (API). 
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PROGRAM SUMMARY 

Program Title: Bill2d 
Journal Reference: 

Catalogue identifier: 

Licensing provisions: GNU General Public License ver¬ 
sion 3 

Programming language: G-l—1-(14) 

Computer: Tested on x86 and x86_64 architectures. 
Operating system: Tested on Linux, and OS X versions 
10.9-10.11. 

Compilers: C-|—1-14 compliant compiler. Tested compil¬ 
ers: Glang 3.7, GGC 4.9-5.2, and Intel G-l—I- compiler 
15. 

RAM: Simulation dependent: kilobytes to gigabytes 
Parallelization: Shared memory parallelization when 
simulating ensembles of systems. 

Keywords: Glassical mechanics; billiards; nonlinear dy¬ 
namics; chaos; transport; diffusion; numerical simula¬ 
tions; molecular dynamics 

Classification: 4.3 Differential equations, 7.8 Structure 
and Lattice dynamics, 7.9 Transport properties, 7.10 
Gollisions in solids, 16.9 Classical methods 
External routines/libraries: 

Boost, CMake, GSL, HDF5; and optionally Google- 
Mock, GoogleTest, and Doxygen 
Nature of problem: 

Numerical propagation of classical two-dimensional sin¬ 
gle and many-body systems, possibly in a magnetic 
field, and calculation of relevant quantities such as 
Poincare sections, survival probabilities, diffusion co¬ 
efficients, etc. 

Solution method: 

Symplectic numerical integration of Hamilton’s equa¬ 
tions of motion in Cartesian coordinates, or solution 
of Newton’s equations of motion if in a magnetic held. 
The program implements several well-established algo¬ 
rithms. 

Restrietions: 

Pointlike particles with equal masses and charges, al¬ 
though the latter restrictions are easy to lift. 

Unusual features: 

Program is efficient, extremely modular and easy to ex¬ 
tend, and allows arbitrary particle-particle interactions. 
Additional comments: 

The source code is also available at https:// 
bitbucket.org/solanpaa/bill2d See README for lo¬ 
cations of user guide, developer manual, and API docs. 
Running time: 

From milliseconds to days, depends on type of simula¬ 
tion. 


1. Introduction 


Hamiltonian systems still offer a wide range of yet 
unexplored territories and interests for the study 
of nonlinear dynamics and chaos. This is evident 
from the recent progress, e.g., in the study of trans¬ 
mission, escape rates, and survival probabilities of 
open systems and recurrences in closed systems P- 
iin], stickiness and marginally unstable periodic or¬ 
bits [Mia, and detailed analysis of the graphene¬ 
like Lorentz gas and related systems dlHll], and 
this list barely scratches the surface. As pointed 
out by Bunimovich and Vela-Arevalo, “Chaos the¬ 
ory is very much alive” [T]. 

Two-dimensional (2D) Hamiltonian systems (i.e., 
with two coordinate dimensions), such as dynami¬ 
cal billiards, have the advantage of being simple, 
but not too simple. That is, they are easy to study 
and visualize, but even a single-particle system can 
display a large variety of complex phenomena of 
Hamiltonian chaos. Many of these systems can 
also be realized as semiconductor nanostructures 
[22 123] • While governed by quantum mechanics, 
these nanostructures have ballistic regions where a 
classical treatment applies to some extent [5] . 

2D billiards (including soft potentials) with 
point-like particles can be simulated either using 
the exact solution, the exact Poincare map, or by 
solving Hamilton’s or Newton’s equations of mo¬ 
tion. The exact solution is often unknown, apart 
from single-particle systems with simple geome¬ 
tries. Therefore, one often resorts to numerical so¬ 
lution of the classical equations of motion. 

In this paper we introduce Bill2d, a code for 
generic, classical 2D systems with a focus on the 
study of chaos and nonlinear dynamics. The code 
can handle single or many particles, hard-wall 
boundaries (as in traditional billiards), particle- 
particle interactions, external potentials, periodic 
systems, a magnetic field, etc. Written in modern 
C-|—k, Bill2d is modular and easy to extend with¬ 
out sacrificing speed or ease-of-use. 

The paper is organized as follows: In Sec. [^we 
describe the class of systems to which Bill2d can 
be applied. In Sec. [^we briefly discuss the numeri¬ 
cal propagation algorithms and in Sec. [^implemen¬ 
tation and structure of the code. In Sec. [^ we give 
a few numerical examples of what can be simulated 
with Bill2d. In Sec. Sec. [Hwe finish with a brief 
summary of the paper. 
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2. Systems 
2.1. Overview 

Bill2d is designed for classical dynamics of in¬ 
teracting particles in two-dimensions. The program 
can deal with both traditional billiards with hard- 
wall boundaries as well as with soft external poten¬ 
tials including also periodic systems. Single- and 
multiparticle systems with a generic form for the 
particle-particle interaction can be simulated, and 
an external magnetic field can be included. Trans¬ 
port calculations are supported by allowing the cre¬ 
ation and removal of particles during the simula¬ 
tion. 

When treating systems without a magnetic field, 
we propagate the particles via Hamilton’s equations 
of motion using symplectic algorithms. Systems 
with a magnetic field are propagated using New¬ 
ton’s equations of motion. 

The program is written using Cartesian coordi¬ 
nates and atomic units, i.e., we set masses, charges, 
and the Coulomb constant to unity. All the par¬ 
ticles are point-like and with equal masses and 
charges (as, e.g., electrons in semiconductors), al¬ 
though the equal mass and charge requirement can 
be lifted with minor modifications to the code. 

To clarify, systems without an external magnetic 
field are described by the fV-particle Hamiltonian 
(7V>1) 

z=l ^ ^ i<j 

where and are the position and momentum of 
the ith particle, H(r) is the total external single¬ 
particle potential, and W{ri,rj) the interparticle 
potential. In principle, potentials V and W can be 
arbitrary, but in the case of attracting singulari¬ 
ties, new propagation algorithms should be imple¬ 
mented. 

When the system includes a magnetic field, we 
use the Newtonian formulation instead of Hamil¬ 
tonian formulation. The equations of motion are 


vector. The first term in Eq. ([^ is the force due to 
the external potential, the second term the inter¬ 
particle interaction, and the third term the Lorentz 
force with the magnetic flux density B perpendic¬ 
ular to the two-dimensional plane. 

The systems can also include hard-wall billiard 
boundaries. The (fully elastic) collision of the ith. 
particle with the billiard boundary is described by 
the transformation 

r, —>■ r, 

(3) 

p, ^ p, - 2 [n(r,) • Pi] n(r,), 

where n(ri) is a unit vector perpendicular to the 
corresponding boundary at position ri. 

2.2. Note on Coulomb interaction and the interac¬ 
tion strength parameter 

When describing Coulomb-interacting systems, 
we introduce the interaction strength parameter a. 
Essentially the Coulomb-potential is given by 

By fixing the total energy and length scale (i.e., 
the size of the billiard table) of the system, we can 
actually study all length and energy scales of all 
geometrically similar systems just by varying a. 

Naturally this trick works only if the system does 
not have external potentials (apart from a hard-wall 
boundary) or a magnetic field. In case of a magnetic 
field and/or external potentials, also those have to 
be adjusted. For detailed derivation of the corre¬ 
sponding scale transformations and demonstration 
of the idea, we refer the reader to Appendix B of 

Ref. [5]. 

3. Algorithms 

3.1. Propagation without a magnetic field 

The propagator, for Hamiltonian systems with 
the Hamiltonian H, is 


( 2 ) 


where ri^^/y is the x or y coordinate of the position 
of the Ah particle and &x/y th® corresponding unit 


B 


dt 


dvi 


dt 


exp(t{-,77}), (5) 

where {•, •} is the Poisson bracket operator. 

The systems we are interested in are described 
by Hamiltonians [Eq. 0 ] that can be split to two 
parts as H = K (p) -|- U{r), where all the momen¬ 
tum dependence is in K, and all the coordinate de¬ 
pendence is in U. As the Hamiltonian is separable. 
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there are plentiful of explicit algorithms for the time 
propagation. This is in contrast to general Hamil¬ 
tonian systems where the symplectic integrators are 
typically implicit and hence computationally more 
demanding. 

The algorithms we have implemented are based 
on split operator schemes, where the propagator, 
exp {At {■, H}) is approximated by some product of 
exp (const • At {•, if}) and exp (const • At {•, Hj). 
The operation of these latter exponentials can be 
carried out exactly - to numerical accuracy - al¬ 
lowing for an easy implementation of split oper¬ 
ator schemes. For example, the velocity Verlet 
scheme [53] can be expressed as 

We have implemented symplectic algorithms of 
orders 2-6 for systems without a magnetic field. 
The algorithms (from Refs. |53H27]) are specified 
in USERGUIDE. 

3.2. Propagation in a magnetic field 

In a magnetic field, the systems are propagated 
via Newton’s equations of motion [Eq. ([^] as we 
are unaware of any efficient symplectic algorithms 
for this class of systems. We have implemented two 
algorithms: (i) the second order Taylor expansion 
scheme developed by Spreiter and Walter [28] and 
(ii) the fourth order algorithm developed by Y. He 
et al. in Ref. [55] . 

The algorithm by Spreiter and Walter has been 
later shown to be an energy conserving split- 
operator scheme [SO]? and the fourth order algo¬ 
rithm is volume preserving (and also preserves the 
total energy extremely well) [55]. These schemes 
explicitly incorporate the effects of the magnetic 
field into the propagation formulas, and allow - in 
principle - an arbitrarily strong magnetic held with¬ 
out the need for a smaller time step. 

4. Implementation 

4 . 1 . Overview 

Bill2d is designed as a modular, object-oriented 
package, and is written in standards-compliant 
C-|—I- following the International Standard 
ISD/IEC 14882:2014(E) Progrcunming Language 
C++ (aka C-l—1-14) [ST]. The package contains 
three binaries: bill2d for single simulation of a 
(many-body) trajectory, bill2d_escape for the 


calculation of escape rates in open billiards, and 
bill2d_diffusion_coefficient for the calcula¬ 
tion of diffusion coefficients in periodic systems. 
All the binaries obey similar inputs; more details 
are given in USERGUIDE. 

The binaries make use of the libbill2d li¬ 
brary, which is created during the compilation pro¬ 
cess. Most of the classes and methods of the 
library are enclosed in the namespace bill2d. 
This library can be used for easy implementa¬ 
tion of additional binaries. A tutorial for devel¬ 
oping new binaries and features can be found in 
docs/developer_manual.pdf, which also describes 
the implementation in detail. 

4 . 2 . Input and initial positions 

The class Parser handles parsing of the input 
either from command line, conhguration hie, or a 
combination of these. The option parsing is imple¬ 
mented with the help of Boost: :program_options 
library. 

The input dehnes the system (number of par¬ 
ticles, magnetic held, potentials, periodicity, etc.) 
and simulation parameters (algorithm, time step, 
simulation time, etc.). The initial positions and ve¬ 
locities of the particles can either be supplied man¬ 
ually, or Bill2d can randomize them. 

The default random initial conditions are calcu¬ 
lated as follows: First the initial positions of the 
particles are randomized either within a given hard- 
wall billiard table, a unit cell, or a user-supplied 
rectangular area. Next, the kinetic energy is dis¬ 
tributed evenly among the particles, and the direc¬ 
tions of the velocities are randomized. Note also 
that implementation of new randomization proce¬ 
dures for the initial conditions is easy. 

4 . 3 . Application programming interface 

The package consists of several classes each de¬ 
signed usually for a single purpose: Table han¬ 
dles collisions with hard-wall boundaries. Datafile 
handles saving of data, etc. The easiest way for a 
developer to create instances of all the necessary 
classes is to use the Parser to get all the input pa¬ 
rameters bundled in an instance of ParsmieterList, 
which can then be passed on to BilliardFactory. 
This factory class can be used to create instances 
of Billiard, which bundles all the other neces¬ 
sary objects together. Furthermore, instances of 
the Billiard class save data automatically during 
simulation and upon destruction. In all simplicity, 
a minimal working example could be as follows. 
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#include "bill2d/parser.hpp" 

#include "bill2d/billiard_factory.hpp" 

using namespace bill2d; 

int main ( int argc , char *argv[] ) { 

Parser parser( argc, argv ); 

auto params = 

parser.getParameterList(); 

BilliardFactory factory( params ); 

auto billiard = factory.get_bill(); 

// Set initial condition 

billiard->prepare (3.0 /* energy */); 

// Get reference to time. 

II This variable is kept up to date, 
const autofe time = 

billiard->get_time (); 

while ( time < 5.0 ) 

billiard->propagate (); 

} // Saving data upon 

// destruction of billiard 

For more complex examples, see files 
src/bill2d_single_ruii.cpp, src/bill2d_ 

escape.cpp, and src/bill2d_diffusion_ 
coef f icient. cpp, and the developer manual 
docs/developer_manual.pdf 

Datafile 

The program saves its data (trajectories, ener¬ 
gies, etc.) to a HDFS-file. The datafile can be ac¬ 
cessed afterwards for analysis with several tools in¬ 
cluding, e.g., the ready made scripts included in the 
Bill2d package, different programming languages 
(at least Python, C, C-|—P, Fortran, and Java) and 
programs such as Matlab and Mathematica. The 
USERGUIDE describes in more detail how to access 
the HDF5-files. 

4-5. Test suite 

The software is bundled with comprehensive unit 
tests, which currently cover most of the program. 
Additional tests can be run with the Python script 
test/propagator_tests .py, which requires the 
user to download reference data from the URL spec¬ 
ified in README. 


5. Numerical examples 

5.1. Trajectories 

A fundamental concept in classical point-particle 
mechanics is a trajectory, i.e., the curve that the 
particle draws in the coordinate space. In Fig. we 
show a single-particle trajectory in the Bunimovich 
mushroom billiards calculated with the bill2d 
binary. The trajectory has been drawn with the 
draw-trajectories script, which is bundled with 
the software package. 

Similarly, one can consider trajectories of a 
many-particle system, as in Fig. which shows 15 
Coulomb-interacting particles inside a circular bil¬ 
liard table. Here several different trajectories share 
the same color. 



Figure 1: Single-particle billiards in a Bunimovich mush¬ 
room. 



Figure 2: Circular billiards with 15 Coulomb-interacting 
particles. 


5.2. Poincare sections 

The phase space can be studied by selecting a 
two-dimensional subspace of the full phase space, 
and visualizing all crossings of a phase space tra¬ 
jectory (or trajectories) with the chosen subspace, 
which is called the Poincare section. 

In Fig. 1^ we show the Poincare section of a 
single-particle elliptic billiards in a magnetic field 
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calculated with the bill2d binary; several initial 
positions are considered in this figure. Here the 
Poincare section shows the collisions with the bil¬ 
liard table, parametrized by the arc length s and 
the tangential velocity ti||. The Poincare section 
consists of several elliptic and parabolic KAM is¬ 
lands (after Kolmogorov, Arnold, and Moser [S31- 
m) that correspond to regular motion, as well as 
a chaotic sea indicated by irregularly distributed 
points. This system has previously been studied in, 
e.g., Ref. m- 



trajectories that can never escape the system are 
not included in the data. 

After a brief initial period, the escape-rate curve 
is found to be a combination of exponential and al¬ 
gebraic functions, which is typical for systems with 
a mixed phase space in the corresponding closed 
system m- In this system the mixed phase space 
(studied, e.g., in Ref. |3S1) is evident from, e.g., the 
regular characteristics of the example trajectory in 
the inset of Fig. 



Figure 3: Poincare section of an elliptic billiard with semi¬ 
major and -minor axes 1 and 3/4 in a magnetic field with 
Larmor radius r = 3. 


5.3. Escape rates in open billiards 

Open billiards have one or many holes in the bil¬ 
liard boundary through which the particle(s) can 
escape the system. The object of interest here is 
the survival probability or the escape rate. The 
survival probability Ps(t) gives the probability that 
the particle has not exited the system before time t, 
and the escape rate Pe{t) gives the probability den¬ 
sity of escape within a narrow time interval around 
t. These two are naturally related by 


dPs{t) 

dt 


Pe{t)- 


(7) 


In Fig. 1^ we show the escape rate for a single 
particle square billiards in a magnetic field, as cal¬ 
culated with the bill2d_escape binary (blue line). 
The square has side length of unity, with a hole of 
length 0.3, and the magnetic field is such that the 
Larmor radius is 0.5075. As a reference, we have 
calculated the escape rate using an analytical map¬ 
ping between boundary collisions (red line), which 
is available for such a simple system. The results 
are nearly identical all the way to very small es¬ 
cape rates (i.e., very rare trajectories). Note that 


Figure 4: Escape rate for a single-particle in square billiards 
in a magnetic field with Larmor radius 0.5075. The result 
demonstrates a combination of algebraic and exponential be¬ 
havior as a function of time. The escape rate computed from 
an ensemble of trajectories using bill2d_escape is shown as 
a blue line, and a reference curve computed with an analyt¬ 
ical collision map is shown as a red line. 


5 . 4 . Diffusion in periodic lattices 

When studying planar billiards in periodic struc¬ 
tures, it is of interest to see how the particle diffuses 
- on the average - in the system. The diffusion co¬ 
efficient is defined as 


D = 


lim 

t—^OO 


{\\r{t)-r{0)r) 

it 


( 8 ) 


where (•) is the average over all possible initial con¬ 
figurations. 

Let us consider a hexagonal antidot lattice 
of scatterers with lattice constant unity (see 
the inset of Fig. [^. Each scatterer pro¬ 
duces a repulsive potential of the form Va{r) = 
(1-I-exp[(||r — RqII —d)/(j])~^. This is defined by 
three parameters: scatterer position R^, radius d, 
and the softness of the potential a. The system 
thus corresponds to a Lorentz gas - an extensively 
studied system in chaos theory - with soft bound¬ 
aries mj. 
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The mean squared displacement (MSD) (|jr(t) — 
r(O)lp) of a particle in the antidot lattice is shown 
in Fig. as a function of time. In this example 
the total energy E = scatterer radius d = 0.15, 
and scatterer softness a = 0.1. The MSD has a 
power law behavior MSD ~ E with 7 « 2, which 
means that the system is superdiffusive, and the 
diffusion coefficient ([^ diverges. This behavior is 
due to large number of trajectories demonstrating 
Levy walk characteristics or even purely ballistic 
behaviour. 

In this example we have simulated an ensemble 
of 5500 trajectories, which results in approx. 2.5 % 
relative standard error of the mean for the MSD at 
t = 4000. The calculation was done with the binary 
bill2d-diffusion_coefficient 

6 . Summary 

The software package Bill2d is a versatile and 
efficient tool to simulate classical two-dimensional 
systems. It supports, e.g., single and many-body 
simulations including varying particle number, in¬ 
terparticle interactions, different external poten¬ 
tials, a magnetic field, periodic boundaries, and 
hard-wall billiard tables. In addition to the ver¬ 
satility, the main strengths of Bill2d are its mod¬ 
ular, object-oriented design, which allows easy im¬ 
plementation of new features, and a clear and ex¬ 
tensively documented source code. 

Bill2d has already been used in several stud¬ 
ies including, e.g., chaotic properties of many- 
body billiards HE], nonlinear dynamics of Wigner 
molecules [35] , and diffusion in periodic lattices [ID] . 



0 1000 2000 3000 4000 

time (a.u.) 


Figure 5: Mean squared displacement for a soft Lorentz 
gas with scatterer radius d = 0.15, softness tr = 0.1, energy 
i? = i behaves as a power law in time (~ t^) with exponent 
7 Ri 2. The inset shows a single trajectory in the system. 


The package can also be readily applied to, e.g., 
classical studies of transport properties in two- 
dimensional structures. 
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